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This paper presents a robust control design methodology for systems 
with probabilistic parametric uncertainty. Control design is carried out by 
solving a reliability-based multi-objective optimization problem where the 
probability of violating design requirements is minimized. Simultaneously, 
failure domains are optimally enlarged to enable global improvements in 
the closed-loop performance. To enable an efficient numerical implementa- 
tion, a hybrid approach for estimating reliability metrics is developed. This 
approach, which integrates deterministic sampling and asymptotic approx- 
imations, greatly reduces the numerical burden associated with complex 
probabilistic computations without compromising the accuracy of the re- 
sults. Examples using output-feedback and full-state feedback with state 
estimation are used to demonstrate the ideas proposed. 

I. Introduction 

Achieving balance between stability and performance in the presence of uncertainties is 
one of the fundamental challenges faced by control engineers. Trade-offs must be made to 
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reach acceptable levels of stability and performance with adequate robustness to parameter 
uncertainty. These trade-offs are explicitly linked to the control engineer’s choice of uncer- 
tainty model as well as how that model is exploited in the synthesis process. Usually, the 
assumed uncertainty model has a profound impact on the performance robustness of the 
closed-loop system. 

Several uncertainty models, such as norm-bounded perturbations, interval analysis, fuzzy 
sets and probabilistic methods 1,2 are typically used. The most commonly used robust con- 
trol methods 3 are //-synthesis and H 00 . In these methods, uncertainty is modeled with 
norm-bounded complex perturbations of arbitrary structure about a nominal plant. This 
treatment is used primarily because it leads to a tractable set of sufficient conditions for 
robust stability, making the approach computationally efficient. These methods are based 
on the most pessimistic value of performance among the possible ones, usually referred to 
as ‘ worst-case 1 . This worst-case performance is usually realized only by a single member 
of the uncertain model set and by a particular input signal. No information is provided 
regarding the likelihood that this worst-case will ever occur in practice. In addition, the 
intrinsic mathematical requirements of the approach usually lead to conservative models of 
uncertainty, over-conservative designs and complicated compensators. 

Probabilistic uncertainty not only defines a set of plants where the actual dynamic system 
is assumed to reside but also associates a weight, i.e., the value of the probability density 
function, to each member of the set. In contrast to conventional robust control methods, 
this ‘additional dimension’ allows the pursuit of robustly optimal solutions in the proba- 
bilistic sense. For instance, reliability-based design searches for solutions that minimize the 
probability of violating design requirements prescribed in terms of inequality constraints. 
Hence, reliability-based control design searches for the compensator that places as much 
probability as possible within the region where the design requirements are satisfied. Notice 
that this allows the search for the compensator with the best robustness for a given control 
structure, e.g. the most robust PID controller, even though the violation of some the design 
requirements for some of the plants in the uncertainty set is unavoidable. 

Synthesis approaches based on random searches 4- ' and stochastic gradient algorithms 8,9 
have been applied to probabilistic robust control. In these studies, random sampling is the 
primary tool for assessing and pursuing acceptable levels of robustness in the control solution. 
On the other hand, asymptotic approximations 10-12 for the estimation of failure probabilities 
have been only used as a control analysis tool. This paper integrates these numerical tools. 
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The main contributions of this paper are as follows: 

1. The use of shapeable failure domains within the reliability formulation. This allows 
the concentration of the random outcome about regions where an improved controlled 
performance is attained. 

2. The formulation of a unified framework where reliability metrics for random variables 
and processes are integrated. 

3. The integrated use of deterministic sampling and asymptotic approximations in a hy- 
brid approach. This approach (i) reduces the computational complexity of the syn- 
thesis algorithm without compromising the accuracy of the results, (ii) eliminates the 
random character of the estimation, and (iii) eliminates the high computational de- 
mands associated with the estimation of small failure probabilities via Monte Carlo 
sampling. These improvements mitigate the high computational demands of existing 
design strategies. 

This paper is organized as follows. Section If presents basic concepts related with control and 
probabilistic uncertainty. Section 111 introduces the reliability metrics for random variables 
and processes to be used throughout the the paper. Realizations to stability, time- and 
frequency-dependent performance metrics are provided therein. Section IV presents the 
hybrid approach used for the estimation of the reliability metrics previously introduced. 
The reliability-based control synthesis procedure is presented in Section V, including robust 
performance considerations and the synthesis algorithm. Two examples are presented in 
Section VI, where a satellite’s attitude control problem and the disturbance rejection in a 
flexible beam are used to demonstrate the method. Finally, some conclusions are stated in 
Section VII. 


II. System Dynamics 

Let p be a vector of random variables used to model the uncertain parameters of the 
system. In this study, p is prescribed a priori by the joint probability density function 
(PDF) / P (p) or equivalently by the cumulative distribution function (CDF) F p (p) a . The set 
of values that p could take, called the support of p, will be denoted as A p . 

a In these expressions, the subscript refers to the symbol used for the random variable while the value in 
parenthesis refers to a particular realization. 
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Consider the probabilistic model J\4( p) of a Linear Time Invariant (LTI) system, where 
the dependence of the model on the uncertain parameters could be non-linear. The reader 
must notice however, that the developments presented herein do not require the system to 
be LTI. The propagation of A p through Af leads to a set of uncertain plant models in which 
the physical system is assumed to reside. The probability of occurrence of a plant within 
this set is fully determined by A4 and p. In a transfer function representation, we will refer 
to G(p) as the uncertain plant and to K{ k) as the compensator, where k is the vector of 
design parameters to be determined. Alternatively, a state space realization of -M(p) leads 
to 


x = A(p)x + B(p)u + F(p)z (1) 

y — C(p)x + D(p)u + E(p)v (2) 

where x is the state, u is the control, z is process noise, y is the system output and v is 
sensor noise. The noise signals are commonly modeled as delta correlated Gaussian white 
noises satisfying E[z] = 0 and E[z(t)z T (t + r)] = S<5(r), where z = [z T ,v T ] T and S is a 
constant spectral density matrix. In what follows, the explicit dependence on p is omitted 
while D is assumed to be zero. 

As a result of uncertainty, important properties used in control design do not hold due 
to the offset between the deterministic mathematical model and the actual dynamic system. 
The effects of parametric uncertainty on the Separation Principle are considered next. For 
the full-state feedback law u = — Gx and a full-order observer with gain L based on the 
expected plant E[Af(p)] (any other deterministic plant such as A4(E[p]) could be used 
instead), the closed-loop dynamics is given by 

x = Ax + Bz (3) 

y = Cx + Ez (4) 


A - BG 

BG 

A — E[A] + (E[B] — B)G 
+L(E[C] - C) 

(B — E[B])G + E[A] — LE[C] 
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where x = [x T , e 1 ] T is the augmented state vector, e = x — x is the estimation error, x is the 
estimation of x, C = [C T |0 r ] T and E = [0 7 IE 7 "] 1 . The vector k is formed by the feedback 
gain G and the observer gain L. Notice that the Separation Principle holds, i.e. A is upper 
triangular, if the deterministic plant used to generate the observer matches exactly the actual 
dynamic system. Uncertainty in the plant makes the Separation Principle unattainable. In 
addition, the random closed-loop poles do not occur at the locations selected for the full-state 
feedback, i.e. poles of the A^i subsystem, nor at the locations for the full-order observer, 
i.e. poles of the A 2j 2 subsystem. 


III. Reliability-Based Metrics 

The propagation of a fixed set of parameters of the plant through conventional control 
analysis tools leads to set of scalar quantities, e.g. closed loop poles, and a set of fields, e.g. 
step responses and Bode plots. The propagation of probabilistic uncertainty through the 
same tools leads to random variables, e.g. random closed-loop poles, and random processes, 
e.g. the step responses become random processes parameterized by time and the Bode plots 
become random processes parameterized by frequency. In this section we first introduce 
reliability metrics for random variables and processes. Such metrics will be used to quantify 
the violation of the design requirements. Specific realizations corresponding to stability, time 
and frequency requirements are then provided. In general, we will use x and x(h) to denote 
a random variable and a random process dependent on p through the plant model. For the 
random process x(h), h refers to an arbitrary parameter such as time or frequency. 

A. Random Variables 

We start by introducing the concept of probability of failure. Let x(p) be the random variable 
of interest and x < x be a design requirement. This event will be referred to as failure. The 
corresponding failure set is given by T = {x \ x G (— oo, a;]}, where the failure envelope x is a 
deterministic quantity prescribed in advance. The admissible domain, namely A = {x \ x E 
(x, oo)}, is the complement of the failure domain. The same type of discrimination can be 
done in the parameter space p by using x(p). The function g(p,x) = x(p) — x, called the 
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limit state function , divides the parameter space in two parts, the domain leading to A, i.e. 
g(p,x) > 0, and the domain leading to J i.e. g(p,x) < 0. Hence, T results from mapping 
the set {p G A p | g(p,x) < 0} through rr(p). In this case, the probability of failure Pf is 
given by 

P f = P[x < x] = f f x [x)dx = f f P {p)dp (5) 

J x<x J g< 0 

Similar expressions can be derived if the design requirement is x > x. A reliability metric 
for x in which constraints from below and above are present is given by 

r x (x, x) = r x (x) + r x (x) (6) 

where 

r x (a ;) = P[x < x\ — F x (x) (7) 

f x (x) = P[x > x\ = 1 - F x (x) (8) 

Notice that r x (x) is equivalent to Equation (5). We will refer to x and x as the envelopes 
of the failure domain T = {x j x G (— oo,a;] U \x, oo)}. Notice that the under-bar and 

the over-bar refer to the bound from below and the bound from above of the admissible 

domain A = {x \ x E (x,x)}. This convention will be used for the reminder of the paper. 
Notice that the mapping of the corresponding limit state function through x(p) leads to the 
failure envelope(s). Hence, there is a direct correspondence between T and g. A sketch with 
relevant information is provided in Figure 1. 

B. Random Processes 

The random process x(h) can be considered as the parameterization of a random variable 
by the deterministic quantity h. The random process x(h) is specified by the set of CDFs 13 
F x (h)( x ,h). For instance, the system output y(t) is prescribed by F y ^(y,t). The evaluation 
of the process at a particular h value, say hi, leads to a random variable prescribed by 
F x (x) = F x (h){x, hi). In general, the support and the percentiles of x(h) depend on h. 
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Figure 1. Sketch on the reliability metric for x. 


In this context, a reliability metric for x(h) is cast as follows 

r,(i,) (s(l>),x(h)) = £*(i,)($(fc)) + r* w (z(fc)) 


( 9 ) 


where 


^2 


^2 


Lxth^yj 1 )) = / P[s(h) < z(/i)]d/i = / F x ( h )(x(h),h)dh 


'hi 

f»/l4 


r»/l4 


r x (h){x{h)) = / P[x(/i) > x(/i)]d/i = / 1 - F x(/l) (a;(h), h)dh 


'h 3 


>h 3 


( 10 ) 

(ii) 


are the costs of violating the constraints ;r(h) < x(/i) and ic(/r) > x(h) respectively. The 
failure envelope functions, namely x(h) and x(h), are deterministic functions that delimit 
the failure domain T = {(x,h) \ (x(h) < x(h) Vh e [hi,/r 2 ]) U (x(/i) > x(/i) Vh G [h 3l h^])}. 
Notice that the admissible domain A is bounded by x(h) from below and by x(h) from above. 
The reader shall realize that Equation (9) is a natural extension of Equation (6). A sketch 
with some of the pertinent metrics is provided in Figure 2. I 11 the top plot, the 1, 25, 75 
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and 99 percentiles' 3 are shown along with the linear failure envelopes x(h) and x(h). In the 
bottom plot, the integrands of Equations (10-11) corresponding to the configuration in the 
top plot are shown. Notice that if the process is contained within the set A the reliability 
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Figure 2. Sketch on the reliability metric for x(h). 


metric r x ^) is zero, meaning that the inequality constraints are satisfied for all parameter 
values in A p . 

C. Realizations 

1. Robust Stability 

A LTI system is robustly stable if all its poles are in the left hand side of the complex plane 
for all possible values of the random parameters. A reliability assessment of stability is given 
by 


P 


U(»M > o) 

. 2—1 


e 


b Recall that the to percentile, given by the x values satisfying (x, h) = to/ 100, defines a line under 
which m% of the probability lies. These lines allow us the visualize the h dependence of the PDF. 
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where s t with i — 1, 2, . . . v is a random pole, 5ft[-] is the real part operator and e is the 
resulting probability of instability. Robust stability is attained if e = 0. Stability can also 
be cast via 


A = max{K[si], K[s 2 ], . . . , 5R[s v ]} (12) 

In terms of A, the probability of instability is given by r,\(0). Robust stability is attained if 
f x ( 0 ) = 0 . 

Several comments are now pertinent. Reaching robust stability may not be feasible for 
the given support A p (even though it is bounded) and the assumed control structure K (k). 
Notice also that the acceptance of a small non-zero probability of instability could be desir- 
able from the performance point of view. For instance, by allowing the right low-probability 
tail of /a (A) to lie on the right hand side of the complex plane significant enhancements 
in the performance of the plants associated with the high probability portions of the PDF 
can be attained. Rather than advocating for the acceptance of the risk that this practice 
implies, we would like to highlight that by allowing small values of e, the trade-off between 
robustness and performance can be quantified. 

2. Time-Domain 

Quite frequently, performance requirements are prescribed in terms of time-domain specifi- 
cations. The propagation of / P (p) leads to random processes for the time responses. Denote 
with x(t) an arbitrary random process with CDF F x u\(x,t). Such process is parameterized 
by time t and the compensator design variable k. The dependence of x(t) on k has been 
omitted for the sake of simplifying the notation. Reliability metrics for relevant processes 
can be cast using Equation (9). For instance, while settling time and overshoot require- 
ments are integrated using r y ^(y{t),y{t)), the control saturation requirement |w| < U leads 
to r ‘ U (t)(—U, U). 

A reliability metric for assessing the effects of noise on the uncertain plant is formulated 
next. The state covariance matrix, defined as Q(t) = E[x(t)x T (t)], is given by the solution 
to the covariance equation 


Q = AQ + QA t + BSB t 


(13) 


subject to Q(0) = Qo- The output covariance, defined as Y (t) = E[y(f)y T (f)], reaches the 
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steady-state Root Mean Square (RMS) value 


y rm s = lim (diag CQ(t)C 

t—> oo V 


1/2 


(14) 


Notice that uncertainty in p makes Yrms a random vector. If Drms is a component of y rms-, 
a reliability metric that penalizes the violation y rms > Vrms i s gi ven by r yRMS (Vrm s ) • 


3. Frequency-Domain 

The propagation of / p ( p) to the frequency domain leads to random processes of the form 
x(lo), fully prescribed by F x ^(x,u>). Here, x{u) is any real frequency dependent metric 
of the feedback loop, e.g. Bode magnitude. This random process is parameterized by the 
frequency oj and the design variable k. A reliability metric for a; (a;) is r x ^(x(u),x(u)). For 
instance, conventional control requirements 14 for disturbance rejection, noise attenuation 
and reference tracking can be cast in terms of the loop transfer function. In terms of the 
loop gain, namely q(cu) = \GK\, r q ^(l) is a metric for low frequency requirements while 
r q (u){q(w)) with q( uj) having a proper roll off, is the metric for high frequency requirements. 


IV. Numerical Estimation 

In general, reliability metrics cannot be evaluated exactly since they involve the eval- 
uation of complicated integrals, usually multi-dimensional, over complex domains. In this 
paper, reliability metrics are estimated using a hybrid approach which combines sampling 
and asymptotic approximations. Such an approach is based on the estimation of failure prob- 
abilities for the random variable. The estimation of failure probabilities for x via sampling 
is given by 




l(xj 6 T) 
n 


(15) 


where X(-) is a binary indicator function that gives one if its argument is true and zero 
otherwise. The subscript in the above expression refers to samples of the dependent random 
variable x. An equivalent expression, where the limit state function is evaluated at the 
sample values of p , can also be used. Usually, Monte-Carlo-Sampling (MCS) 4 ~ 6,15 is used to 
generate the required samples. 
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A. Hammersley-Sequence-Sampling (HSS) 

HSS generates representative deterministic samples of / P (p). The error of approximating 
an integral by a finite sample of the integrand, e.g. Equation (15) instead of Equation 
(5), depends on the uniformity of the points used to generate the samples rather than on 
their randomness. This has motivated the development of deterministic sampling techniques 
such as HSS, where the distribution of points is optimized. The n Hammersley samples 
are generated by transforming the n Hammersley points irq with % — 1, 2, , . . n through the 
inverse CDF of the uncertain parameter 

P* = X ( m i) (16) 

The Hammersley points can be generated 16 easily. HSS requires far fewer samples 17 than 
conventional MCS for a given confidence level. Improvements in the convergence rate of the 
estimated first two order moments by a factor of three to one hundred 18 have been reported. 
In addition, if HSS is used to generate the samples for Equation (15) the estimated value 
of the failure probability is deterministic. In contrast, MCS leads to a random value for 
Pf unless an infinite number of samples is used. This is especially noticeable if n is small. 
The random character of the estimation can only be mitigated by increasing the number 
of samples, which incidentally increases the computational demands of algorithms based on 
MCS. Therefore, HSS not only leads to more accurate estimations than MCS for a given 
number of samples but also eliminates the random character of the results. 

Figure 3 shows a comparison between HSS and MCS. In the top, n = 200 points on the 
unit hypercube are shown. In the bottom, the corresponding samples for / p ( p) = f a (a)fb(b ), 
where f a (a ) = 1V(0, 1) and fb(b ) = B( 3,2) with At, = [0,1] are displayed. Here, N and B 
denote a Gaussian and a Beta distribution. Substantial differences in the uniformity of the 
points and in the clustering of the samples are observed. 

B. First-Order-Reliability-Method (FORM) 

FORM 11 uses an asymptotic approximation for the estimation of failure probabilities. In 
the process, p is transformed into the standard normal uncorrelated space q. If p = T(q) = 
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Figure 3. Points and samples via MCS and HSS. 


F p 1 (F q (q)), Equation (5) is equivalent to 

Pf= I / q (q)dq 

Jg{ T(q))<0 

FORM approximates the domain g (T(q)) < 0, by a half-space fitted to the true domain at 
the point of maximum probability density. This approximation leads to 

Pf ~ $Hlq*ll) (17) 

where q*, called the Most Probable Point (MPP), is given by the solution to the constrained 
optimization problem q* = argi||q||, p(T(q*)) = 0. In this expression, $ refers to the CDF 
of a standard normal random variable. Notice that the rotational symmetry of / q (q) leads 
to the one-dimensional approximation in Equation (17). The MPP does not exist when the 
probability of failure is zero or one since the equality constraint p(T(q*)) = 0 cannot be 
satisfied. Even though FORM is extensively used in structural engineering, its application 
to controls has been limited to stability 10 analysis. The use of FORM to estimate failure 
probabilities related to A might lead to non-smooth limit state functions. This occurs when 
crossings between any pair of limit state functions corresponding to the system poles take 
place. 
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C. Hybrid Approach 

Sampling based techniques can readily be used to estimate probabilities of failure using 
Equation (15). However, high computational demands in the evaluation of ay = x(pj) can 
preclude their practicality especially when Pf ~ 0. Examples of this can be easily found 6 . 
On the other hand, methods based on asymptotic approximations, such as FORM, provide 
good approximations when Pf is small. This is clear since for failure probabilities away 
from zero and one, the slow decrease in F q (q) near the MPP and the geometrical difference 
between the true limit state function and its linear approximation contribute a bigger error 
to the FORM approximation. 

In this paper, a hybrid approach which combines HSS and FORM is used to estimate 
probabilities of failure. In order to identify the numerical tool that best suits the task at 
hand, a coarse and computationally-efficient estimation of Pf is first generated using HSS. 
Snch estimation is then compared with a reference, namely the reference failure probability 
p, a user-defined value set in advance. The comparison between the coarse estimate and p is 
used to determine if either FORM or HSS are used to generate the more accurate estimation. 
Assume that two sets of Hammersley samples of / p ( p) are available. One set has n \ samples 
and the other one has n 2 samples, where ri 2 n 1 . For a given failure domain T and a given 
reference failure probability p, proceed as follows 

1. Estimate the Pf using Equation(15) and the set of ri\ samples. 

2. Recalculate Pf as follows. If the estimated value is greater than p use Equation (15) 
with the set of n 2 samples. If the estimated value is less than p use FORM. 

The refinement performed in Step 2 might not always be necessary. Situations in which 
this is the case are provided next. Since reliability metrics for random processes are heavily 
dependent on the larger values of the probabilities of failure that compose them, (see the bot- 
tom plot of Figure 2), refining the estimation of the relatively small values is inconsequential. 
Furthermore, if the reliability metrics are used to calculate the cost function of an optimiza- 
tion problem, more accurate estimations are not needed when the assessment resulting from 
using the coarse estimate of Step 1 denotes a poor control design, e.g. ^(0) 0. 

c Wang et al. 5, 15 make the approach computationally viable by using a single random variable to model 28 
uncertain parameters. Same authors 6 require 25000 samples to determine a sufficiently small 95% confidence 
interval. 
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The above procedure applies to the random variable x. Extensions to random processes 
are easily attained. For the random process x(h), generate e samples in the h domain as 
follows 

% — 1 

hi zrihmax hmin) “ 1 “ hmim ^ 1 , 2 ... 6 

e — 1 

where h max and h m i n are the bounds of the interval of interest, e.g. h max = ^2 and h nun = h\ 
in Equation (10). Reliability metrics for the e random variables x(hi), are estimated via the 
hybrid approach and then used to form the integrands in Equations (10-11). For the sake of 
clarity, we will refer to the e samples in the h domain as the e partitions. 

V. Control Synthesis 


A. Problem Formulation 

The formulation of the control design problem from a reliability perspective is as follows. 
For a given set of design requirements, plant model, compensator structure and uncertainty 
model we would like to find the compensator parameters for which the resulting probability 
of violating the controlled system requirements is minimized. Such requirements combine 
stability and performance specifications in time and frequency domains. Notice that this 
statement refers to the excursion of the outcomes into the failure domains. Performance 
improvements regarding the outcomes within the admissible domains will also be considered. 
Such improvements are attained by dynamically shaping the failure domains during the 
minimization of the reliability cost metrics. This topic is considered next. 

B. Performance Improvements 

The reliability metrics in Equations (6-9) are applied using a fixed failure set T . In this 
form, a reliability analysis cannot assess the system’s performance in the regions where the 
design requirements are satisfied, i.e. the union of the admissible domains A associated 
with all the design requirements. Since the portion of the random outcome lying on the 
admissible domain A might end up being substantially larger than the portion lying on 
the failure domain T . a reliability-based approach with fixed failure envelopes do not have 
control over the bulk portion of the PDF, which is the portion that dictates the most likely 
overall performance. 
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The ideas behind the approach to be proposed will be introduced with an example. Let 
x(k) be the stationary RMS value of an error signal. Usually, we would like to find k such 
that x is as close as possible to zero. Uncertainty in the plant makes x a random variable. 
Let x be the failure envelope associated with a design requirement, i.e. T = {x \ x E [cc, oo)}. 
The minimization of f x {x) leads to a reliability optimal compensator. Suppose there exist 
multiple designs leading to r x (x) = 0. These designs however differ in how well the resulting 
PDF of x spreads over the admissible domain A = {a; j x G [0,3;)}. The concentration of 
f x (x) about zero is an indicator of the overall performance. Say, k x leads to r x (x/ 2) = 0 and 
k 2 leads to r x (x) = 0. Since none of these two designs violate the design requirement x > x, 
a reliability analysis cannot establish that the compensator with parameters k x has a better 
global performance than the one which uses k 2 . 

By minimizing the reliability metrics and simultaneously enlarging the failure domain(s), 
the whole random variable/process can be concentrated about regions where an improved 
system performance is achieved. This is attained by parameterizing both the failure envelopes 
and a failure size penalizing function with an additional design variable. This variable will 
be denoted as c. For the RMS example above, the minimization of J = r x (c ) + c, where the 
design variable is d = [k, c] and c G [0, x\. leads to solutions which integrate reliability and 
performance considerations on a single formulation. Notice that the value of J for k x is less 
than the one for k 2 if c G \x/2,x). 

In general, we will refer to the augmented reliability metric as the sum of a reliabil- 
ity metric from Section III and a penalizing term. Augmented reliability metrics for the 
random variable x and the random process x{h ) take the form r x (x( c),x(c)) + 7 X ( c ) and 
r x (h.)(x(h, c ),x(h, c))+7 x(h) (c) respectively. The penalizing functions 7 x (c) and 'y x (h.)(c) must 
be proportional to the size of the admissible domain A. In addition, they must be built such 
that the minimization of the augmented metric does not lead to inadmissible solutions, e.g. 
r x = 1 and 7 x = 0. If r < e is desired, use a monotonically increasing function satisfying 
7 G [0, e]. For the RMS example above, the minimization of the augmented metric leads to 
f x < e if 7 x(c) = ec/x for c G [0, x]. 

C. Synthesis Procedure 

A step-by-step procedure to reliability-based control synthesis is presented next. 

1. Determine the plant model and the control structure. First principles and classical 
deterministic approaches to compensator design can be used. Identify the set of pa- 
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rameters that have a strong impact on the plant model. Use sensitivity information 
and engineering judgment to select the set of uncertain parameters p. At this stage, 
the parametric plant model, e.g. G(p), and the control structure, e.g./l(k), must be 
fully determined. 

2. Generate the probabilistic parameter model / p ( p). Use engineering judgment and 
experimental data if available. 

3. Determine the number of HSS samples n \ and n 2 to be used. Follow the guidelines 
provided in Section IV-C. Use Equation (16) to generate the sample sets of / P (p) for 
both ni and n 2 . 

4. Cast the violation of the design requirements in terms of reliability metrics as in Equa- 
tions (6-9). Recall that specific realizations for stability, time and frequency require- 
ments were provided in Section III-C. Use these metrics to compose the reliability cost 
vector r. This step requires determining the failure domain T corresponding to each 
component of r. 

5. Determine which failure domains will remain fixed during synthesis and which ones 
will be dynamically shaped. Let d be the design variable. When the failure domains 
are fixed d = k. When failure domains are shapeable d = [k 7 , c T ] T . 

6. Build a penalizing function y(c) for each of the components of r whose failure domain 
is not fixed. Follow the guidelines in Section V-B. Update the components of the 
reliability cost vector r by adding the penalizing functions and parameterizing the 
failure envelopes. 

7. Solve the multi- objective optimization problem 

J = min{r r Nw} (18) 

d 

where N is a diagonal normalization matrix such that the components of r 7 N are 
between zero and one, and w is composed of non-negative weights. The matrix N is 
used to facilitate the weighting of the components of r. 

Each cost function evaluation used in the search for the optimal reliability-based design 
d* requires a reliability analysis. This analysis is done by calculating the reliability 
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metrics contained in r using the hybrid approach. This task requires forming the 
closed-loop Equations (4-14) and performing typical control studies such as finding 
closed loop poles, time responses and Bode plots. 

During optimization, the following procedure is suggested in order to focus most of the 
computational effort toward the assessment of better designs. First, calculate the cost 
function using ri\ samples for e\ partitions. This implies that only the first step of the 
hybrid approach is applied to all reliability metrics. This Erst assessment, denoted as 
Ai, should be computationally efficient. If A\ shows that d is a good design relative 
to the ones already evaluated by the optimizer, perform the refined assessment A 2 . 
The assessment A 2 is carried out by using e 2 partitions and a adjustable value for the 
reference failure probability p. The adjustment of p is done to prevent inconsequential 
calculations as it was mentioned in Section 1V-C. If the particular failure probability 
is to be estimated via HSS use n 2 samples and Equation (15). If the particular failure 
probability is to be estimated via FORM use Equation (17). If Ai indicates that d is 
not a good design, A 2 is not carried out. This two-fold analysis is applied to all designs 
the optimizer evaluates in the search for d*. 

Implementing the dual assessment described above avoids the inconsequential refinement of 
the reliability metrics. 

D. Optimization and Reliability 

Due to the nature of the reliability metrics in r, the cost function J(p,d) might not only 
have plateaus, i.e. , there could exist a design d and a non- zero perturbation 5 such that 
J(p,d) = J(p,d + 5), but might also have a discontinuous gradient. 

The use of sampling in the estimation of probabilities makes the cost function discon- 
tinuous at every point of the design space. Let J(p,d) be an estimation of the actual cost 
J(p,d). For any design d and regardless of the number of samples, there always exists a 
perturbation 5 such that J(p,d) = J(p,d + <5). This situation is aggravated, i.e. bigger 
perturbations can be found, when a smaller number of samples is used or when Pf is close 
to zero or one. 

The discontinuous nature of the estimated value of J must be taken into account when 
selecting a numerical optimization method to solve Equation (18). In the examples to follow, 
the resulting non-convex non-continuous optimization problem is first solved using Genetic- 
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Algorithms (GA) for a fixed number of generations. Since GA is based on a random search, 
the twofold procedure described above is particularly convenient. After the fixed number of 
generations is reached the GA solution is refined using the Nelder Mead Simplex algorithm, 
which is a local non-gradient based search method. 

VI. Numerical Examples 

The synthesis procedure of Section V-C is applied herein. A textbook satellite attitude 
control problem is considered first. Then, disturbance rejection for a flexible beam is pre- 
sented. If p G M m , the parameters used for A\ are ri\ = 75 m and e\ = 90. For A 2 we use 
n 2 = 500 m and e 2 = 180. For the sake of comparison, the examples present the solution 
to deterministic versions of the problems for which the expected value of p is used. Such 
problems and the corresponding solutions are referred to as the nominal ones. 

A. Attitude Control 

Accurate satellite pointing in the presence of large thermal gradients and mass losses for 
uncertain initial conditions is desired. A simple rotational model of two bodies connected 
with a flexible boom leads to 


AA + b(9\ — 0 2 ) + k{6\ — 0 2 ) — u 
A A + b(0 2 — A) + k(9 2 — A ) = 0 


where A and 9 2 are the deflection angles, A and J 2 are moments of inertia, k is the equivalent 
stiffness, b = ayJkJlA is the equivalent damping coefficient and u is the applied torque. The 
variable a is used to model the changes in damping caused by thermal variations. We 
assume that A = 0.1 since mass losses only affect A- The non-collocated sensor-actuator 
pair resulting from using y = 9 2 leads to the SISO transfer function for the plant 


G(p) 


k + bs 

A.'A'S 4 + 6(</i + J 2 )s 3 + (A + J 2 )ks 2 


(19) 


Variations in the operating conditions and the ignorance on the initial conditions are modeled 
using p = [A, e, k, 9, 9] T , where 9 = A(0) = A(0). The following output-feedback control 
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structure is assumed 


( 20 ) 


. , k\ + k 2 s + k 3 s 2 + k 4 s 3 
k 5 + k e s + k 7 s 2 + /c 8 s 3 

The joint PDF that describes the uncertainty in p is given by the independent random 
variables listed in Table 1, where U and B refer to Uniform and Beta distributions. Notice 
that the Beta distribution has four independent parameters, two of them are the conven- 
tional arguments and the other two are the support bounds. Performance requirements 


Table 1. Uncertainty model. 


Jl 

A j! = [0.8, 1] 

fji{Ji) = U{ 0.8,1) 

a 

A a = [0.03,0.2] 

f a (a)=B( 0.3, 0.2) 

k 

A fc = [0.09, 0.4] 

f k (k) = B( 5,5) 

9 

A 0 = [— vr/2,7r/2] 

fe(9) = B( 5.2, 5.2) 

9 

A<j = (—15, 15) 

f e {9) = B{ 2.5, 2.5) 


on the system’s closed-loop stability, settling time, over-shoot, and control usage for a 
step response and on the magnitude of the loop transfer function lead to r = [r\(A), 

*?(«)(?M)] T » where = \ GK \ is the lo °P § ain - 
The failure envelopes to be used are A = 0; y(t) = *—1.25 'Hit) + 2.2H(t — 70) for t E [0, 80], 
where H is the Heaviside function; y(t) = 1.2577(f) — 0.2 H(t — 70) for t E [0, 80]; u(t) = —0.5 
for f E [0,25]; u(t) = 0.5 for f E [0,25]; q(u) = 0.75/a ; for u E [10~ 6 , 0.2] and q(u) = 1 for 
uj E [1, 10 2 ]. These failure domains lead to the normalization matrix N = diag{ [1, 80, 25, 0.2— 
10 -6 , 10 2 — 1]}. Note that overshoot and settling time constrains are integrated using y(t). 
Fixed failure domains are assumed. This leads to d = [k\, k 2 , k 3 , k 4 , k 5 , k 6 , k 7l k 8 ] T ■ The solu- 
tion of the optimization problem in Equation (18) requires calculating r for multiple control 
designs. Recall that for each design, the hybrid approach of Section V-C is used to generate 
the coarse assessment A 4 and eventually the refined assessment A 2 . For A 4 , use ri\ = 375 
samples, ei = 90 partitions and Equation (15). If A 2 is required, use e 2 = 180 partitions, 
the reference failure probability p = 0.01 and the results of A\. For the failure probabilities 
to be estimated via HSS use n 2 = 2500 samples. In this example, FORM is used to calculate 
fR(0) for p = 0.01. Before presenting the results from the above procedure, the deterministic 
problem based on the expected value of the parameter is considered. 
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1. Nominal Compensator 


A baseline compensator for the nominal plant is designed by standard pole placement tech- 
niques such that large stability margins are attained. This practice results in the nomi- 
nal compensator parameters k = 10 6 [0. 0108, —0.3271, 0.1192, 0.0092,1-8835,2.1305,2.2276, 
0.9308] t . A reliability analysis of the nominal compensator using the probabilistic uncer- 
tainty prescribed by / P (p) leads to r T N = [0,0.2485,0.227,7.08 x 10~ 5 ,0]. This vector 
indicates that the closed-loop system is robustly stable, i.e. Fa(0) = 0, but the time re- 
sponses are unsatisfactory. The CDF of A as well as the time evolutions of the output and 
the control signals are shown in Figures 4-6. The sudden variation in the slope of the CDF of 
Figure 4 is the result of a change in the closed-loop pole that determines A. The considerable 
disparity between A(E [p] ) and E[A(p)] shows that the nominal problem is not a meaningful 
representative of the probabilistic behavior. Figures 5 and 6 show the time evolution of the 
random signals by indicating the 1, 10, 20, 30, 40, 50, 60, 70, 80, 90 and 99 percentiles. 
In Figures 5 and 6, the percentiles, the nominal fields and the failure envelopes are shown. 
Dotted lines are used to indicate the the failure envelopes. It is interesting to see how the 
PDFs expand, e.g. Figure 6 at 2.5 and 8 seconds, and contract, e.g. Figure 6 at 4 and 
16 seconds, in a oscillatory manner. This information can be used to determine the time 
periods when the effects of uncertainty are more noticeable. 


KE[pl) Efl.(p)l 



Figure 4. A for the nominal compensator. 


20 of 30 


°v 

5 1' 

0 —• 

5 

0 / \ 
I#' ■ 


\ \\\ 
' V \ ' 


\ 99 percentile 




Percentiles 
y(t) at E[p] 




■/■;yr 

•a:.:::-'"- 

1 1,1 1 percentile 

till 

i i * i i 




Failure envelopes 



Figure 5. y(t) for the nominal compensator. A zoom is shown below. 



Figure 6. u(t) for the nominal compensator. 


2. Reliability -based compensator 

The synthesis procedure listed above leads to d* = 10 6 [0.0405, 0.1267, 0.2422, 0.0320, 0.5244, 
1.0057, 1.2263, 0.6560] r and r T N = [3.13 x 10“ 4 , 0.0521, 0.0918, 1.33 x 10~ 4 , 0] for which the 
weighting vector w = [500, 1, 1, 1, 1] T was used. A probabilistic analysis of this compensator 
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leads to Figures 7-10. Notice that despite of the increased variability of the dominant closed- 
loop poles resulting from this compensator (Aa in Figure 7 is about three times larger than 
the one shown in Figure 4), stability is compromised with only 3.13 x 10~ 4 probability. 
Notice that a substantial improvement in the performance is achieved by trading-off a very 
small margin of the probability of instability. This improvement can be seen after comparing 
Figures 5-6 with 8-9. Better robust stability characteristics could be attained by increasing 
the weight in w for the corresponding component of r. Recall that reaching zero probability 
of instability might be unfeasible. From Figure 9 we see that for all possible parameter 
values and initial conditions the process u(t) stays between the ±0.5 range with more than 
0.8 probability after 6 seconds. Figure 10 shows that uncertainty mostly affects the damping 
and the location of the resonant frequency. Violations of the low frequency requirement are 
completely avoided. Overall, the performance resulting from d* is substantially better than 
the one resulting from k. 

During optimization, 157 random variables were used to evaluate c for the coarse assess- 
ment A\ . Such task takes 23.6 seconds when performed on a Pentium 111 1795MHz with 
512MB of RAM. Notice that the CPU time associated with A 2 depends on the initial con- 
ditions used to find the MPPs. For this assessment, HSS was used for 628 random variables 
and the hybrid approach was used for robust stability. This task took 102 seconds. 


E[X(p)] X(E[p]) 



1 

Figure 7. A for the reliability-based optimal compensator. 
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Figure 8. y(t) for the reliability-based optimal compensator. A zoom is shown below. 



Figure 9. u(t) for the reliability-based optimal compensator. 


B. Disturbance Rejection for a Flexible Beam 

The second example will focus on a reliability-based disturbance rejection solution for a 
flexible beam test article with both physical and modal parameter uncertainties. The system, 
displayed in Figure 11, consists of a very flexible thin aluminum blade, approximately one- 
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Frequency [rad/s] 


Figure 10. Bode plot of the loop gain for the reliability-based optimal compensator. 


meter long, attached at its base to a hub motor. The hub motor is the control actuator for 
the system. At the tip of the beam, there is a reaction wheel that serves as a disturbance 
generator. The test article has nine sensors that may be used in any combination for either 
feedback or performance output monitoring. The finite element method is used to model 
this system by utilizing Euler-Bernoulli planar beam elements. A complete description of 
the flexible beam test article 19 is available. 

For this paper we study a SISO model of the system in which the input u is the hub motor 
torque and the measured output y is the tip velocity. The tip reaction wheel disturbance 
is modeled by passing a Gaussian white noise process through a second-order linear low- 
pass filter, with parameters £ = 0.8 and u) n = 2007T rad/s. The first five modes of the 
clastic structure are used to build a state space realization of the plant. This, in addition 
to the disturbance model leads to a open-loop system where x G M 12 , u G M. and y G 1. 
The uncertain parameters are the Young’s Modulus E (Pa), the density p (Kg/m 3 ) and the 
damping ratios of the retained vibration modes. This set leads to p = [E, p, £i, £ 2 , £3, £4, £s] T , 
whose components are assumed independent. The corresponding PDFs are given in Table 
2. The mean value of the parameters E[p] is set to coincide with the parameters in the 
finite element model. These mean values were chosen to match experimental data, while 
the supports of the distributions were set according to reasonable ranges of variation. The 
shapes of the PDFs were arbitrarily set. Performance requirements on stability and the 
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Figure 11. Flexible beam test article. 


Table 2. Uncertainty model. 


E 

A e = 10 10 [5.226, 7.839] 

h(E)=B( 5,5) 

P 

A p = [2280, 3420] 

fM = B( 3,3) 

6 

A Sl = [0.08,0.12] 

> 

To 

to 

6 

A 6 = [0.0252, 0.0378] 

Mb) = B( 2,2) 

6 

A& = [0.02, 0.03] 

> 

CO 

CO 

JO 

to 


A ?4 = [0.0304, 0.0456] 

MU) = B( 2,2) 

£5 

A 6 = [0.02, 0.03] 

= B( 2,2) 


output RMS lead to r = [r\(0),r y RMs(yRMs)] T ■ Two control design problems are considered 
in this example, output feedback and full-state feedback with a full-order observer. 

1. Output Feedback 

First, a third-order compensator with the same structure of Equation (20) is considered. A 
baseline compensator is designed such that the tip velocity RMS for the nominal plant is 
minimized. This results in i/rms = 0.025 m/s and k = [—0.0052, —0.2589, —16.1462, —5.007, 
0.0004, 0.0009, 0.2519, 0.0598] T . A reliability analysis of the nominal compensator for the 
uncertainty model in Table 2 leads to r A (0) = 0.068 and r y KMs(0-05) = 0.010. 
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For the reliability based-design, a shapeable failure domain for the RMS component of r 
is assumed. This leads to the cost vector r = [r A (0), r y RMs(c) + 7 y RMs] T , where c G [0,0.05] 
and 7 v rms — c. The corresponding design variable is d = [k T , c] T and the normaliza- 
tion matrix is N = diag{ [1, 1.05]}. The weighting vector w = [20, 1] T leads to d* = 
[-0.0072,-0.3506,-22.2811, 0.0145,0.0006,0.0019,0.3761,2.8 x 10" 8 , 0.0305] T , f A (0) = 0 
and r. y /jA^s(0.0305) = 4.33 x 10 -3 . A probabilistic analysis of d* leads to Figures 12-13. Fig- 
ure 12 shows that the the whole random variable yRMS is moved towards zero, by virtue of 
the non-fixed failure envelope. Figure 13 shows Bode magnitude plots of the disturbance to 
output transfer function, namely T zy , Notice that differences in the low-frequency portion of 
the diagram have a bigger impact on the RMS value. In this Figure, Bode magnitude plots 
for the nominal compensator and the reliability optimal compensator are shown. The 1, 10, 
20, 30, 40, 50, 60, 70, 80, 90 and 99 percentiles and the nominal fields are superimposed. 



Figure 12. Tip velocity RMS for the output feedback solution. 


2. Full-State Feedback and Full-Order Estimation 

Second, full-state feedback and a full-order observer define the structure of the compensator is 
considered. Hence, the feedback gain G, the observer gain L and the RMS failure envelope c 
are the design variables. Recall that the separation principle does not hold due to uncertainty 
in the plant. The resulting closed-loop dynamics is given by Equations (3) and (4). Notice 
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Figure 13. Bode diagrams of T zy for output feedback. 


that although the observer is deterministic, all the closed-loop poles are random. As before, 
a baseline compensator for the nominal plant is first designed by minimizing the RMS value, 
which results in Drms = 0.011 m/s. Using y RM g = 0.05 m/s, this compensator leads to 
Fa(0) = 0.232 and r yR Ms( 0.05) = 0.002. 

The reliability-based synthesis approach for the same setup used in the output feedback 
example leads to r A (0) = 0 and ^^^5(0.0139) = 3.6 x 10~ 3 . Due to the large number of 
elements in d*, only c = 0.0139 m/s is provided. The probabilistic analysis of d* leads to 
Figures 14-15. Compared with Figure 12, Figure 14 shows considerably more variability in 
the Bode magnitude plot as well as a significant reduction in the damping of the first mode. 
Since there is no conservatism in the selection of the nominal plant, i.e. G(E[p]) is not the 
most difficult plant to control, the optimal deterministic values for the nominal plant don’t 
have to bound the resulting supports, e.g. values in A yRM g may be less than y R MS — 0.011 
m/s. 

ft is interesting to notice that even though d* leads to a robustly stable closed-loop 
system in Equation (3), the full-state feedback subsystem A M and the full-order observer 
subsystem A 2j 2 have a non-zero probability of instability. This indicates that the Separation 
Principle artificially reduces the design space. While robust stability was achieved in both 
output feedback and full-state feedback solutions, the latter led to a better performance. 
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y RMS 

Figure 14. Tip velocity RMS for the full-state feedback solution. 



Figure 15. Bode diagrams of T zy for full-state feedback for full-order observer. 


VII. Conclusions 

This paper proposes a reliability-based control synthesis method for systems with proba- 
bilistic uncertainty. Control synthesis is performed by solving a multi-objective optimization 
problem where the probability of violating stability and performance requirements is min- 
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imized while the failure domains are simultaneously enlarged. Including dynamically sha- 
peable failure domains leads to improvements in the global controlled system performance 
that could not be pursued by a reliability formulation with fixed failure domains. In addi- 
tion, the integrated use of asymptotic approximations and deterministic sampling in a hybrid 
approach proved to considerably relax the high computational demands of the synthesis al- 
gorithm. Examples of attitude control of a simple satellite model and disturbance rejection 
of a flexible beam are used to demonstrate the method. 

References 

1 Crespo, L. G., “Optimal performance, robustness and reliability based designs of systems with struc- 
tured uncertainty,” Proceedings of American Control Conference, Vol. 5, Denver, CO USA, June 2003, pp. 
4219-4224. 

2 Laughlin, D. L., Jordan, K. G., and Morari, M., “Internal model control and process uncertainty - 
mapping uncertainty regions for SISO controller design,” International Journal of Control, Vol. 44, No. 6, 
December 1986, pp. 1675-1698. 

3 Zhou, K. and Doyle, J. C., Essentials of Robust Control, Prentice Hall, Upper saddle, New Jersey, 

1998. 

4 Marrison, C. and Stengel, R., “Design of Robust Control Systems for a Hypersonic Aircraft,” Journal 
of Guidance, Control, and Dynamics , Vol. 21, January-February 1998, pp. 58-63. 

5 Wang, Q. and Stengel, R. F., “Robust Nonlinear Control of a Hypersonic Aircraft,” Journal of Guid- 
ance, Control, and Dynamics, Vol. 23, No. 4, July-August 2000, pp. 577-585. 

6 Wang, Q. and Stengel, R. F., “Robust control of nonlinear systems with parametric uncertainty,” 
Automatica, Vol. 38, 2002, pp. 1591-1599. 

7 Crespo, L. G. and Kenny, S. P., “Reliability-based control design for uncertain systems,” AIAA Journal 
of Guidance, Control, and Dynamics, Vol. 28, No. 4, 2005. 

8 Polyak, B. and Tempo, R., “Probabilistic robust design with linear quadratic regulators,” Systems and 
Control Letters, Vol. 43, 2001, pp. 343-353. 

9 Lagoa, C. M., Li, X., and Sznaier, M., “On the Design of Robust Controllers for Arbitrary Uncertainty 
Structures,” IEEE Transactions on Automatic Control, Vol. 48, No. 11, November 2003, pp. 2061-2065. 

10 Spencer, B. F., Sain, M. K., Won, C.-H., Kaspari, D. C., and Sain, P. M., “Reliability-based measures 
of structural control robustness,” Structural Safety, Vol. 15, 1994, pp. 111-129. 

11 Rackwitz, R., “Reliability analysis, a review and some perspectives,” Structural Safety, Vol. 23, 2001, 
pp. 365-395. 

12 May, B. S., Probabilistic Robust Control: Theory and Applications, Ph.D. Dissertation, California 
Institute of Technology, Pasadena, CA USA, 1997. 


29 of 30 



13 Crespo, L. G., “Probabilistic formulations to robust optimal control,” 45th AIAA Structures, Structural 
Dynamics and Materials Conference, Palm Springs, CA USA, April 2005, pp. 1-21, AIAA Paper No. 2004- 
1667. 

14 Skogestad, S. and Postlethwaite, I., Multivariable feedback control, John Wiley and Sons, Chichester, 
England, 1996. 

15 Wang, Q. and Stengel, R. F., “Searching for Robust Minimal-Order Compensators,” Journal of Dy- 
namic Systems, Measurement and Control, Vol. 123, June 2001, pp. 233-236. 

16 Hokayem, P., Abdallah, C., and Dorato, P., “Quasi-Monte Carlo Methods in Robust Control Design,” 
IEEE Conference on Decision and Control, Maui, HA USA, December 2003, pp. 2435-2440. 

17 Kalagnanam, J. R. and Diwekar, U. M., “An Efficient Sampling Technique for Off-line Quality Con- 
trol,” Technometrics, Vol. 39, No. 3, August 1997, pp. 308-319. 

18 Diwekar, U. M. and Kalagnanam, J. R., “Efficient Sampling Technique for Optimization under Uncer- 
tainty,” American Institute of Chemical Engineering Journal, Vol. 43, No. 2, February 1997, pp. 440-447. 

19 Kenny, S. P., Optimal rejection of nonstationary narrowband disturbances for flexible systems, Ph.D. 
thesis, Massachusetts Institute of Technology, Cambridge, MA USA, February, 2002. 


30 of 30 



